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Universality of generalized bunching and efficient assessment of Boson Sampling 
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It is found that identical bosons (fermions) show generalized bunching (antibunching) property 
in linear networks: The absolute maximum (minimum) of probability that all N input particles are 
detected in a subset of K. output modes of any nontrivial linear M-mode network is attained only by 
completely indistinguishable bosons (fermions). For fermions 1C is arbitrary, for bosons it is either 
(i) arbitrary for only classically correlated bosons or {ii) satisfies K. > N (or K, = 1) for arbitrary 
input states of N particles. The generalized bunching allows to certify in a polynomial in N number 
of runs that a physical device realizing Boson Sampling with an arbitrary network operates in the 
regime of full quantum coherence compatible only with completely indistinguishable bosons. The 
protocol needs only polynomial classical computations for the standard BosonSampling, whereas an 
analytic formula is available for the Scattershot version. 

PACS numbers: 42.50.St, 03.67.Ac, 42.50.Ar 


Introduction. - Optical networks with photons have be¬ 
come a growing research field with potential applica¬ 
tion in quantum computing [H-Q. Boson Sampling (BS) 
idea 0, a non-universal but near-future feasible device 
aiming at the Extended Church-Turing thesis (ECT), 
followed by spectacular experiments with optical net¬ 
works of growing size and number of photons [iSl: 
can be a way for benchmark demonstration of quantum 
supremacy. A BS device with a random but known M- 
mode network, TV ~ 30 single photons at known input 
modes for M iV^ would achieve this ultimate goal 
‘M- However, the very computational complexity of BS 
111 fl^ requires an exponential in N number of runs of 
such a device and computation of an exponential num¬ 
ber (at least 0{N^)) of classically hard permanents to 
prove BS by comparing an output distribution with the¬ 
oretical probabilities. Quantum supremacy demonstra¬ 
tion thus faces a big challenge of maintaining the TVth 
order quantum coherence in a BS device with iV ^ 30 for 
an exponential number of runs. Though a distribution 
claimed to simulate BS, as the uniform one [T^, can be 
exposed with only polynomial number of runs of a device 
[U (see also Ref. Q), finding such a protocol for a given 
sampler, e.g., the one with distinguishable particles, is a 
hard open problem. 

The proof of BS being exponential both in the number 
of runs and computations does not prevent efficient veri¬ 
fication of the very source of quantum supremacy of BS, 
i.e., the full iVth-order quantum coherence in a device 
with N bosons, where unwanted distinguishability, pho¬ 
ton losses, and higher photon numbers are the leading 
adversary factors in optical setups nm. Such an as¬ 
sessment of BS may require only a polynomial number 
of runs and polynomial classical computations. Zero- 
transmission laws in the Fourier network [l^ and sta¬ 
tistical benchmarking probed this path, but an as¬ 
sessment protocol applicable to an arbitrary network is 
an open problem. The test of Ref. [l^ does verify the 
TVth order coherence but only in a single Fourier network 
(not posing a threat to the ECT), whereas the statisti¬ 
cal method of Ref. [l^ can only assess the second-order 


coherence. Already with relatively small networks exper¬ 
imentalists have to resort to either Bayesian methods or 
circumstantial evidence 1M3- 


(i) (II) 



FIG. 1: (Color online) Assessment protocol for a BS device. 
Step (I): the source is checked for output with only one parti¬ 
cle per mode. Step (II), the statistics of all N input particles 
to land in K, output modes of a M-mode network U is gathered 
by on-off type detectors in complementary L = M — 1C modes: 
The maximum probability is attained only by completely in¬ 
distinguishable bosons (and is close to 1 for NL M). 


Generalized bunching and efficient assessment of BS.- 
We call a protocol that uses only polynomial classical 
computations to certify in a polynomial number of runs 
that a BS device operates at the full quantum iVth-order 
coherence compatible only with completely indistinguish¬ 
able bosons an efficient assessment protocol for BS. Effi¬ 
ciency requires an independent network certification: An 
efficient assessment of BS not using certification of a net¬ 
work matrix in a regime different from BS allows for loop¬ 
holes, due to a combined effect of network and distin¬ 
guishability errors (see the Appendixb This conclusion 
agrees with the results of Refs. [H, [ij]. An assessment 
protocol for a BS device with an arbitrary network can 
only be based on a network-independent (universal) abso¬ 
lute maximum or minimum of some probabilit y at tained 
only by completely indistinguishable bosons [46|. The 
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discovered generalized boson bunching is such a prop¬ 
erty: The absolute maximum of probability of detecting 
all N input particles in K. output modes of an arbitrary 
nontrivial M-mode linear network is attained only by 
completely indistinguishable bosons (i) over arbitrary in¬ 
put states of particles for /C > (and /C = 1) or {ii) for 
arbitrary 1 < /C < M — 1 over only classically-correlated 
bosons. The minimum probability is attained by com¬ 
pletely indistinguishable fermions, with that of distin¬ 
guishable particles lying in the middle between the two. 
The generalized bunching/antibunching is found by ex¬ 
ploiting the discovered equivalence between probability of 
detecting all input particles in a subset of output modes 
of a linear network and an eigenvalue problem for positive 
semi-definite (p.s.d.) Hermitian matrix (Eq. Q below). 

The generalized bunching can be viewed as an N- 
particle generalization of the Hong-Ou-Madel (HOM) ef¬ 
fect [13 to arbitrary M-mode networks. Known gen¬ 
eralizations of the HOM effect to many-particle multi- 
mode setups are reported for special Bell multiport net¬ 
works |18l - [^ . whereas bunching to a single mode of a 
random network [2l| has probability decaying as ^ e~^ 
for M > N'^ [ 2 ^. A similar effect is boson tendency to 
form clouds in output modes in some networks Q. In 
contrast, the generalized hunching can be observed in a 
polynomial in N number of experimental runs in an arbi¬ 
trary (nontrivial) quantum network. Though in the dilute 
limit, M > N'^, the generalized bunching/antibunching 
effect wanes as —>■ 00, it nevertheless allows for an ef¬ 
ficient assessment protocol of a BS device, sketched in 
Fig. [U requiring, besides only a polynomial number of 
runs of a device, polynomial classical computations (see 
below). Moreover, analytical results are derived for the 
Scattershot version of BS [1^ . 

Description of partially distinguishable identical parti¬ 
cles in a linear network.- We consider an Af-mode lin¬ 
ear quantum network with N single identical particles 
in an arbitrary internal state at input modes fci,..., k^ 
(ordered products of operators are assumed in case of 
fermions) 


the key role (23 - [^ . e.g., particles of one species with 
piint) an^isynimetric under permutations emulate behav¬ 
ior of the other species . The probability formula 

of an output configuration m = (mi,..., tum) [H, [ 2 ^ 
(applicable also to fermions, see the Appendix) reads 


p(m) 



N 

t,(T^Sn ol—1 


_ ( 3 ) 

where Ii,...,In are output modes with multiplicities 
(mi,...,mM), and a complex-valued function J(cr) of 
a € Sn is defined as 


J(a) = e(a)Tr(p<“>P.,), £(.) = { 

( 4 ) 

where Pa'[\a=i'^\ja) = na=i® li<T-i(a)) is an operator 
representation of a in the Hilbert space . 

Identical particles are called completely indistinguish¬ 
able if is symmetric under permutations, = 

s(a) (e.g., particles in the same internal state), whereas 
particles with orthogonal internal states are distinguish¬ 
able, (see also Refs. [ 2 ^, |3l|). The 

trace of in the symmetric subspace with 

Sn = {l/N\)Y.aP-^ i-e-, d{J) = Tr{^ivpi“‘i} = 
M 'TIiu^Sn is a suitable measure of the partial 

indistinguishability 1^ and also in other quantum 
information applications of identical particles [33|. 

Note that J{a) of Eq. ([4]) is a p.s.d. function of cr S 
Sn, i.e., for any complex-valued function z[a) we have 
^ z*{ai)J{<7ia2^)z{<72) > 0, while J{I) = 1 for the 

CTl ,(72 

identity permutation I. Therefore, there is a factorizing 
function 0{a) such that (see the Appendix) 


E d*{T)9{Ta), ^ |6»(cr)p = l. (5) 

t^Sn (7^Sn 


N 

p = ^g.|vE-.)(vE-.|, n“l.,j0), (1) 

i j a=l 

where > 0, 3 = (ji, ■ • • ,/Ar), Ej = 

1, and a mode operator aj, j (and below 6^^) creates a 
particle in an input (output) mode k and an internal basis 
state |j) £ H (e.g., a basis function of spectral shape 
of a photon). A unitary network with matrix U, Fig. 
mb), relates the input Okj and output bij modes: ^ = 

'l2fLiUk,ibjj. The internal state of identical particles, 
defined as 

N 

i j a=l 

governs their behavior in a linear network. Symmetry 
properties of under the symmetric group Sn play 


Importantly, any normalized p.s.d. function J(cr) cor¬ 
responds to an input state of single particles (see the 
Appendix), i.e., the whole set of input states with N sin¬ 
gle particles can be equivalently represented by the whole 
convex set of normalized p.s.d. functions Sn —>■ C. 

Probability to detect all particles in K, output modes.- 
The probability for all N input particles to gather in 
/C output modes, say 1,...,/C as in Fig. mb), reads 
Pn{J) = where m = (mi,..., 0 ..., 0). 

Defining an A^-dimensional p.s.d. Hermitian matrix H, 
built from the submatrix of U on the rows fci,...,fcjv 
and columns 1,...,/C, and the corresponding (A^!)- 
dimensional Schur power matrix n(M) indexed by ele¬ 
ments of iSjv, 

K N 

n.,r(i^)= n^-(«)Ma), (6) 

1^1 
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we obtain pn{,J) in the form (see the Appendix) 

Pn{J) = Y. = Y ^*Wnr,r'0(T'). (7) 

(J^Sn t.t'^Sn 

Therefore, Pn{J) is a convex combination of eigenvalues 
of n(77). For the completely indistinguishable particles 
0 h'i)((j) = e(cr)/-\/]V!, whereas for distinguishable parti¬ 
cles = Saj- Eq. (O gives for these ideal cases: 

N 

= Per(i7), p^P = det(i7), ( 8 ) 

Q!=l 


for bosons, fermions, and classical particles, respectively. 
Well-known inequalities (see the Appendix) result in the 
following order 


(F) , (d) , (B) 

Pn <Pn <Pn ■ 


(9) 


Eq. dH) seem to suggest that pP and pP are the ab¬ 
solute maximum and minimum of pn{J)- To establish 
such a property one has to show that they are the unique 
maximum (minimum) eigenvalues of n(iJ). A famous 
result of Schur [s^ states that the smallest eigenvalue 
of n(iJ) is det(i7), hence, the generalized antibunching 
is an universal attribute of the completely indistinguish¬ 
able fermions (a unique minimum for K, > N). On the 
other hand, the maximum eigenvalue of n(Ff) is gener¬ 
ally unknown. The permanent-on-top conjecture (POT), 
stating that universally it is per(i?)j_proven for A^ < 3 
has turned out false for > 5 (3^ . 

Thus, conditions on input state and/or network are 
needed to ensure the maximum probability being at¬ 
tained only by the completely indistinguishable bosons. 
Eor /C = 1 Eq. © gives pn{J) = d{J)N\llY 
with the maximum for an arbitrary network at d{J) = 1 
or J(ct) = 1 (see also Refs. [U, [H, [s^). Numerical 
simulations [47j| with random p.s.d. Hermitian matrices 

/ D\ 

reveal that p)^ ' = per(iJ) is not the maximum probabil¬ 
ity only ior N > 5 particles in 2<IC<N—1 output 
modes (such a state of A^ > 5 bosons is necessarily a 
state of non-classically correlated particles, see below). 
More importantly, when }C > N, a unique eigenvector 
SP{a) = l/vVI corresponds to per (if), i.e., the maxi¬ 
mum probability of N particles to be detected in 1C > N 
output modes is attained only by the completely indis¬ 
tinguishable bosons. Thus, bosons show the generalized 
bunching property in }C > N output modes. 

For only classically correlated bosons, a unique max¬ 
imum of pn{J) is attained only by the completely in¬ 
distinguishable ones for all 1 < /C < M — 1. In¬ 
deed, a classically-correlated internal state can be ex¬ 
pressed as a convex combination of pure states, i.e., 
p{int) ^ 0 ... 0 for some ar¬ 
bitrary states \(j)j^) £ H and Vj > 0, = 1- The 

corresponding J-function reads 

( 10 ) 

j a=l 


Setting we get from Eqs. dl]), ©, and 


PNpp=Y^iP^<H-Gp, (11) 

j 

where the dot stands for the Hadamard (by-element) 
product. Thus, the permanental version of Oppenheim’s 
inequality [s^, stating that for two p.s.d. Hermitian ma¬ 
trices H and G (for Ga,a = 1), per(i/ • G) < pei{H) 
would imply the claimed result in this case. Using ma¬ 
trices H that violate the POT conjecture, it was checked 
that per(iJ • G) < per(iJ) for any N random states 
|(/)i),..., |^i>Ar), with at least two linearly independent. 

Assessment protocol of a BS device.- The aver¬ 
age {pn{J)) over Haar-random networks gives an 
idea of quantitative features of the generalized bunch¬ 
ing/antibunching effect. By the unitary invariance of the 
Haar measure, the average probability (p)^ ' ) depends 

only on the ratio of considered output configurations (see 
the Appendix): 

/C(/C±l)-...-(/C±jVTl) 

' M{M ± 1) •... • (M ± AT T 1) ’ 

here (and below) the upper (lower) signs stand for bosons 
(fermions). For distinguishable particles there is no exact 
result, but for M 1 it can be shown that (see the 
Appendix) 


(*>» > = (^ 


N r 


1 + 0 




KM 


The average probability ratio becomes 


{pPP 


1 + 0 


1^1 


n N-1 

n 

-I i=l 


l±l/K 

l±l/M' 


(13) 


(14) 


For NL M, where L = M — K, the detection prob¬ 
ability is close to 1: Pp'^'^) = ~ L/{M ± 1)] = 

1 — 0{LN/M), whereas the r.h.s. of Eq. (03 gives 
Pp^'^)/Pp) ~ 1±LA^(A^— 1)/(2M^). In this case one 
needs R ^ M'^ / runs for the ratio (fT4ll to surpass 

the statistical error 0(l/-\/i?) in experimental data. But 
the ratio in Eq. (03 is attained only by the completely 
indistinguishable bosons (fermions), hence, it is a reli¬ 
able witness, detectable in polynomial number of runs, of 
their complete indistinguishability during propagation, 
i.e., that no decoherence process has contributed to out¬ 
put statistics. Therefore, we have an efficient protocol 
for assessment of a BS device with an arbitrary network, 
Fig. [TJ The only known protocol for a BS device with an 
arbitrary network 14 1 , experimentally verified dis¬ 
criminates the BS and uniform distributions. Our pro¬ 
tocol discriminates against all other than BS samplers, 
which are physically realizable with particles in a linear 
network, including the classical {Ma), the fermion {J-a), 
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the random-classical ( Ba) samplers of Ref. 0 , and the 
random-phase bosons [l^, [ig. 

The assessment protocol has two stages. At stage 

(I) , Fig. [ijl), by using photon-number resolving de¬ 
tectors (e.g., by cascading bucket detectors [s^) one 
checks that sources produce N single photons. At stage 

(II) , Fig.[T](II), employing the universality of generalized 
bunching, one verifies experimental statistics against the 
probability ([8]) using L = M — K. bucket detectors. 

The protocol requires only one matrix permanent 
= per(iJ) of H in Eq. ([5]) (a single set of /C modes is 
used) to an error e = 0{N~'^) for some k > 0 (statistical 
error in experimental data for a polynomial number of 
runs). For iV ^ 1 in the dilute limit M = 0(A^^+'^) 
with (5 > 0 and L = 0{N) it can be shown (using 
that we select K. modes arbitrarily) that only polynomial 
in N computations Cn are required (see the Appendix; 
in this case (pjvf^) = 1 — 0{N~^) and (pj^^ — Ptv^) = 
One can estimate that, on average over all 
choices of K. output modes, Cn = 0{N ^) (e.g., setting 
K = 2[1 -|- 5] allows one to distinguish the quantum and 
classical cases). 

The protocol applies also to the Scattershot BS [^ . 
recently experimentally tested 0 , which uses heralded 
single photons in N random input modes in each run. 
The probability describing experimental statistics of a 
Scattershot BS is well approximated by (p^^) (fl^ al¬ 
ready for a few hundred runs of such a device (see Fig. 

1 181 below), i.e., no computations required. 

Stage (I) is designed to expose all attempts to bypass 
the universality property using inputs with variable num¬ 
ber of particles per mode (not required under certified in¬ 
put). Indeed, an input with any distribution of the com¬ 
pletely indistinguishable bosons between M input modes 

p = '^Pn\n){n\, Pn>0, ^p„ = l, (15) 


also the sampler Ba of Ref. 0, a “mockup distribu¬ 
tion of BS” physically realized by distributing N uncor¬ 
related distinguishable particles randomly over N input 
modes, with the probability of such a particle to land 
in an output mode I being p{l) = The 

probability of single particles at input is N\/N^ ^ e~^. 


N 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

L 

2 

2 

3 

4 

5 

5 

6 

7 

7 

8 

9 

9 

10 

11 

11 

12 

13 

14 

M 

5 

8 

13 

18 

25 

32 

41 

50 

61 

72 

85 

98 

113 

128 

145 

162 

181 

200 


TABLE I: Network size M and L = M — 1C as functions of N. 
Here K, is selected by maximizing the ratio of Eq. Gl under 
the condition that (p^^) > 0.25 (note that tC > N). 



FIG. 2: (Color online) The analytical quantum average prob¬ 
ability Eq. (I12II (dots on the solid line) and the ap¬ 
proximation Eq. m of the classical one (dots on the 

dashed line) vs. numerical averaging with 1000 Haar-random 
networks (circles). The squares give the quantum and clas¬ 
sical probability for the Scattershot BS (estimated with 500 
runs) with a randomly chosen network for each value of N. 


where ni um = N, has the Haar-average proba- 

bility equal to ) (see the Appendix). For example, 
take the random-phase bosons of Refs. 00 , i.e., an 
input where each boson is in a coherent superposition of 
S input modes (and in an internal state |(())) described by 
operator A{9) = S~^ with random phases 

01,...,% 00. The density matrix of this input reads 


c- 27r 

_ jS-iy.S^ ^ 

{S + N-l)l^J^J 2 tt 
J-r 0 


[At(0)]^|O)(O|[A(%" 


{S + N-1)1 




(16) 


where n = (rik^, ■ ■ ■ , n^s ) > + ... + nks = N. A source 

of Ps with S = N is exposed at stage (I) by a vanishing 
probability of an input with one particle per occupied 
mode, for A' 1 scaling as ^ 4“^. Stage (I) exposes 


Fig. [TH] gives numerical results, where M = [^N^] 
(integer part), and L = M — Kis obtained by maximizing 
the ratio (fTT)) for > 0.25 (see Table H]). 

The distinguishability error 1 — d[J) Ri (1 — F){N — 1) 
[ 3 ^ in a BS device (F is the mean fidelity of indistin- 
guishability of photons) can be assessed from an exper¬ 
iment using the first-order approximations (see the Ap¬ 
pendix): 


Pw ^ -Pn{J) 


1 - d(J) 
A- 1 


Np^P 


jL 

dx 


per{H{x)}^^i 


(pT) - (PNiJ)) « [1 - d{J)]^{p^^l,), (17) 

where Ha,p{x) = Sa,pHa,px+{l-Sa,p)Ha,p (per{A(a;)} 
is a polynomial in x of order A). The second law is valid 
for M A^ and applies to the Scattershot BS. 

Difference between an experimental and the theoretical 
probability would reflect presence of other 
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errors in a device. How the BS regime is affected by er¬ 
rors in a network matrix can be estimated beforehand. 
The aim of such a certification, besides eliminating the 
possibility of loopholes (see the Appendix), is to guar¬ 
antee that network errors would be below an acceptable 
level in the BS regime. The theoretical basis is provided 
in Refs. Besides an unwanted distinguishability 

and network matrix errors, errors reported in BS experi¬ 
ments include higher-order photon numbers, esti¬ 

mated at stage (I) of our protocol, and the photon losses, 
which can be directly estimated at a network output (see 
also Ref. HI). A BS device with a fixed-ratio of lost 
photons is believed to be hard to simulate on a classical 
computer UM with some progress in proof Since 
a lossy linear M-mode network is equivalent to an 2M- 
mode unitary one, with one half of output modes being 
unaccessible “loss channels” (see the Appendix), our as¬ 
sessment protocol applies also to linear lossy networks. 
In this case in Eq. (jH]) the proper non-unitary network 
matrix U must be used. It can be experimentally char¬ 
acterized with only classical coherence (45l| . 


Conclusion.- We have discovered the generalized bo¬ 
son bunching and fermion antibunching in linear net¬ 
works and proposed an assessment protocol for BS veri¬ 
fying in a polynomial number of experimental runs that 
a BS device with a random linear network operates at 
the full TVth-order quantum coherence compatible only 
with N completely indistinguishable bosons, i.e., the very 
physical origin of its quantum supremacy. The proto¬ 
col requires only polynomial classical computations for 
the standard version of BS, whereas for the Scattershot 
version (with better prospects for scalability) analyti¬ 
cal results are available. In general terms, the gener¬ 
alized bunching is a generalization of the famous HOM 
effect, revealing the complete indistinguishability of N 
photons in an arbitrary (nontrivial) M-mode network, 
which may find other applications whenever linear quan¬ 
tum networks are used. 
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Appendix 


I An efficient assessment protocol for BS on 
an uncertified network matrix allows for loop¬ 
holes 

Verification of Boson Sampling (BS) with a random 
network without knowledge of a network matrix is an 
ill-defined problem [U, In view of this, a natural ques¬ 
tion arises: Is it somehow possible to verify the full iVth 
order quantum coherence compatible only with N com¬ 
pletely indistinguishable bosons and a (special) network 
matrix simultaneously in a polynomial number of runs in 
the same BS regime on such a device? This may seem 
indeed possible due to presence of a high symmetry in a 
network matrix, leading to very distinct features in out¬ 
put distribution. Let us call such an assessment of BS 
the holistic assessment. 

However, careful analysis reveals loopholes in an effi¬ 
cient assessment protocol if the latter is used as a holistic 
one. Let us analyze in detail an assessment protocol for 
BS based on the Fourier network Q (a highly symmet¬ 
ric Bell multiport), where a large fraction of the output 
configurations are forbidden ]w symmetry for completely 
indistinguishable bosons Q. Let us recall the neces¬ 
sary details of the assessment protocol proposed in Ref. 
Q. One considers Fourier network with M = N'p modes 
(where p > 2), i.e., given by the matrix 

For single particles in the input modes 1 < ki,..., k]\f < 
M with the following cyclic symmetry 

+ NP-^ (19) 

where, obviously, 1 < fci < Np~^, only those output 
configurations 1 < li,... Jn < M which satisfy 

N 

Y^la^qN, qe Af, ( 20 ) 

a=l 

i.e., have the sum of mode indices divisible by N, are re¬ 
alized with completely indistinguishable bosons in such 
a network Q. Moreover, the number of realized output 
configurations is only a fraction 1/N of the total num¬ 
ber of all a priori possible ones, which allows efficient 
verification of the full fVth order indistinguishability by 
counting the number of violations of the forbidden events, 
on the order of 1 — 1 /N of the total number of events Q . 

The above described protocol is subject to loopholes, 
if it is understood as a holistic one, i.e. when no cer¬ 
tification of a network matrix in a different regime is 
performed [l^. Loopholes appear due to a combined 
effect of network and distinguishability errors. Let us il¬ 
lustrate this by exposing a loophole in the following case: 


N = 2Ni and p = 2 (for simplicity). Denote Mi = Nf. 
Since no information on a network is assumed, instead of 
(|18l) one cannot rule out the following one (P is a 
permutation) 


U = P 


p(Mi) 0 p(Mi) pt 


p{Mi) 0 p{Mi) pt 

( 21 ) 

where [/ is a block-structured matrix with two diagonal 
blocks, each containing two Mi-dimensional Fourier ma¬ 
trices (fTSl) whose rows and columns are permuted by P 
in such a way that each odd (even) row and column in 
the two blocks of U of size 2Mi corresponds to one and 
the same submatrix 

Consider now N bosons divided into two groups of 
Ni completely indistinguishable bosons, whereas bosons 
from different groups being distinguishable (e.g., bosons 
in group i = 1,2 are in an internal state {(fi) and 
(^i|02) = 0). Assume that such bosons are launched 
into U cm observing the cyclic symmetry dm), where 
the bosons in the first Ni input modes of U are from the 
same group (i.e., completely indistinguishable). Since N 
is even, the parity of an input mode ka of U is the same 
as that of ki. Thus, the first (second) group of Ni com¬ 
pletely indistinguishable bosons are launched in one of 
the two submatrices of the respective block, in a 

similar cyclic symmetry as in Eq. jni) for its proper 
mode indices and N substituted by Ni. By Refs. [3,13 
and the block-structure of U, the allowed output config¬ 
urations must satisfy conditions similar to Eq. (EOl) now 

(i) (i) 

for the proper output indices, say, 1 < Z) ,..., < Mi, 

of that particular submatrix F^^^\ i.e.. 


Ni 

Y, = q^'^Ni, G A?, i = l,2. (22) 

Of—1 

Let us now verify that for a cyclic input (I19|l all al¬ 
lowed output configurations in a network U cm, with the 
above described partially distinguishable bosons, belong 
to the allowed set A of the Eourier network F^^'> (flSl) 
with completely indistinguishable bosons (where there is 
an exponential number of such |A| ^ First, let us 

assume that ki is even. The indices 1 < li,... ,1m < M 
of output modes of U, where bosons from one of the two 
groups can end up, are derived from the output indices in 
the corresponding submatrix F^^^\ satisfying Eq. (1221) . 
by the following two rules (respectively, for the first and 
second blocks of U) 

L = 2l^^\ Im,+<, = 2Zi2) + 2 Mi, a = l,...,Ni. (23) 

Since N = 2Ni Eq. (l20ll is satisfied by such output mode 
indices 

N 

= (g(l) +g(2) +Mi)fV. (24) 
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Similarly, when fci is odd we have a similar relation for 
the output mode indices 

= 21^^^ - 1, In,+c. = 2/i") - 1 + 2Mi, a = 1 ,... ,iVi, 

(25) 

with a similar conclusion. Thus, all of the exponentially 
many allowed output configurations in U belong to A. 
Though we have considered an even number of bosons, 
a similar example can be devised for an odd number as 
well. 

The considered combination of a matrix U and par¬ 
tially distinguishable bosons is just a special simple case 
of many such combinations of a network and partially dis¬ 
tinguishable bosons (see also Ref. @) with the allowed 
output configurations being an exponential subset of A. 
For example, network U dSD relates its output modes 
to input modes of the same parity, but a slightly modi¬ 
fied network consisting of U cm followed by a “noisy” 
network implementing a random shift of all modes by 
1 (mod M) with probability 1/2 will have an exponen¬ 
tial number of allowed output configurations all belong¬ 
ing to A and no parity symmetry Due to many 

such possibilities, there is no way to guess, not know¬ 
ing a network matrix U ^ what kind of partic¬ 

ular feature distinguishes an exponential in N subset of 
A which corresponds to such U and partially distinguish¬ 
able bosons. An experimentalist assuming that a network 
matrix is close to dni) and input bosons close to 

being completely indistinguishable, attempting a holis¬ 
tic assessment by the above protocol, runs only the BS 
regime. After a polynomial number of runs of a device, 
which in reality has one of the alternative possibilities 
of network matrix U ^ F^’^'> and partially distinguish¬ 
able bosons, the experimentalist would conclude that the 
input is fully compatible with the completely indistin¬ 
guishable bosons and the matrix with F^^\ since with 
only a polynomial number of runs one cannot tell from 
an experimental statistics that only some exponential in 
N part of A is actually realized. 

The conclusion is that, to close all possible loopholes 
in an efficient assessment protocol of BS, it is absolutely 
necessary to independently certify a network matrix, e.g., 
using an efficient method of Ref. @ requiring only clas¬ 
sical coherence. Besides eliminating the possibility of 
loopholes, such a certification would also guarantee that 
network errors are below an established acceptable level 
when BS is run on such a network. The theoretical ba¬ 
sis for a network assessment is provided in Refs. [3, Q; 
where the effect of network matrix errors on output prob¬ 
ability distribution of a BS device with such a network is 
estimated. 


II Probability of N particles to gather in K. out¬ 
put modes of a linear network 


bosons/fermions at input of a linear network U (for 
bosons this result was derived in Refs. i i). We 
will consider simultaneously both species (in case of 
fermions an order of creation and annihilation opera¬ 
tors is assumed). A more general input is assumed, with 
0 < n-fe < particles per input mode k (fermions have 
linearly independent internal states in each input mode). 
A general input state of configuration n = (ni,..., um) 
reads 

P(n) = 


with gi > 0, % = 1) and 




\/Ai(n) 


N 

(0 TT .t 






JO), 


(27) 


where a basis state |j) G TL in the internal space is in¬ 
troduced (e.g., a basis function of spectral shape of a 
photon), j = and /r(n) = Hfeli Per¬ 

mutation symmetry (anti-symmetry) of creation opera¬ 
tors for bosons (fermions) allows to chose expansion co¬ 
efficients symmetric (anti-symmetric) with respect 
to the Young subgroup Sn = 0 ... <8) Sum of the 

symmetric group Sn, where corresponds to permu¬ 
tations of internal states of particles in input mode k 
between themselves. Such coefficients are normalized by 

The probability of an out put configuration m = 
(mi,..., ttim) is given as [1, II^ 

p(m|n) = Tr(p(n)T>(m)), (28) 

where p is the input state Eq. (1^ and T>(m) is the 
detection operator [3, [l^ (|0) is Fock vacuum state) 


27(m) = 


M(m) 


E 


N 




.Q!=l 


r N 


n 


(29) 


One can evaluate the trace in Eq. (1^ by first express¬ 
ing the input mode operators in Eq. (E5)) through the 
output ones using and then employ 

the following identity (see also Refs. EilI3) 


( 0 | 


N 


n <>■..) 


a=l 


N 




|o) 


N 


= E n' 

ct^Sn a—1 




(30) 


Substituting Eq. (IMl) and (l2^ into Eq. (l28)) and using 
Eq. (1501) in the two inner products one obtains the prob¬ 
ability of an output configuration m in a linear network 
U in the form 


p(m|n) 


1 

/r(m)^(n) 


N 


r,a^SN a—l 


Let us first recall principal steps in derivation of 
the output probability distribution of N identical 


( 31 ) 














where are output modes, 1 < < M, with 

multiplicities (mi,... ,mM), whereas function J{cr) and 
the internal state are 


i 


and 


(32) 

j «=i 


J(cr) = e(cr)Tr(p(“‘^P„), e{a) 


1, Bosons, 
sgn((T), Fermions, 


(33) 

where PaUa=i®\ja) = na=i®li<T-i(a)) is the operator 
representation of a in . Note that the permuta¬ 
tion symmetry (anti-symmetry) of the input state (1261) - 
(ITZl) for bosons (fermions), i.e., = pi“‘^Pn- = 

e(7r)pi*”‘i for any tt G 5n, implies that 


with mic+i = ... uiM = 0. We have 
Pn{J) = ^':p(m|n) = ^ ^ ^ 


Zl = l lr\T — l 


N 


Xn) 




T,(7GSi\r 


N 


n(n) X! "^(^ ) n ^a,fT'(a) - , . X! 

^ cr'eSw a=l ^ <y&SK 


(36) 


P( 


- ^ 0*{T)e{r')Ur.y. 

T,T'(iSN 


where we have transformed the sum over output config¬ 
urations m into that over output mode indices li,... ,li\j 
with the combinatorial coefficient defined the 

following matrices 


J{(J'k) = J{-K(j) = J{cr), Vtt G Sn. (34) 

A network input consists of completely indistinguish¬ 
able bosons (fermions) if the corresponding J-function 
reads X‘^)(u) = £((t) @,13 (in case of fermions p.(n) = 
1). This case allows one to completely neglect the in¬ 
ternal degrees of freedom. Probabilities at a network 
output are expressed through the usual matrix perma¬ 
nent and determinant, respectively. The simplest case 
of completely indistinguishable particles consists of all 
particles being in the same internal state |(^), giving 

p(™*) = (|(/.)(<^|)®'^. 

The other limiting case, which may be identified as 
the classical case, since the output probabilities are the 
same as in the case of classical particles, corresponds to 
a “block-structured” J-function (see also Refs. 10) 

(35) 


Function J{a) of Eq. (1^ appears when the inter¬ 
nal states of identical particles from different input 
modes become orthogonal: Tr{p(™*)Po.} = 0 for a ^ 
iSn, whereas (by the symmetry of Cj) we always have 
e((7)Tr{p(*”‘)Po.} = 1 for (7 G (i.e., distinguishable 
particles from the same input mode cannot be discrim¬ 
inated by a linear network from the completely indis¬ 
tinguishable bosons). Note that the subgroup 5n acts 
as identity on the indices ki,... ,kN of matrix U in Eq. 
m, thus the sum over 5n in Eq. (1551) cancels ^(n) in 
the denominator in Eq. m, resulting in the familiar for¬ 
mula for the probability in the classical case, expressed 
through the matrix permanent of doubly stochastic ma¬ 
trix with elements \Uki\'^. Obviously, for single particles 
at input {uk < 1) we have = S^j. 

The total probability of detecting all N input particles 
at a preselected (and fixed) set of 1 < /C < M output 
modes, say, the first K, modes, is a sum of p(m|n) (1511) 


N 


^ Xr(X = n 


<y{oi) ,t(Q:) ’ 


(37) 






reordered the product nLl X(a),.(o) 

Y[a=i ^a,Ta-^{a)^ and defined a' = Ta~^ (the sum 
over T cancels the factor 1/A^!). 

Note that for a p.s.d. Hermitian H, as a principal 
submatrix of the tensor product matrix the Schur 

power matrix is a p.s.d. Hermitian as well. 

Consider now the limit cases. In case of the com¬ 
pletely indistinguishable particles (since 

j(*<i)((j) = e((T)). Eq. ((551) gives for bosons and fermions: 


(B) 

Pn 


per(Fr) 

M(n) 


(F) 

Pn 


det(i7)J^^n),i5 


(38) 


where per(A) = 0^=1 ^ct,CT(a)- For classical particles 

(e.g., distinguishable bosons or fermions from different 
input modes) we get from Eq. (1551) that 


0(^)(cr) 


1 


^ ^ ^{7,775 
7rC«Sn 


(39) 


which results in 


N 

pf= n (40) 

a=l 


The following order of the limit-case probabilities can 
be easily established 


(41) 


valid for arbitrary number of particles per input mode. 
Indeed, for a p.s.d. Hermitian matrix H, which we rear¬ 
range in a block-matrix form 


/ ^(1.1) ^(1.2) \ 
77(2.1) ^(2.2) j ■ 


(42) 









9 


the following inequality is known for the matrix determi¬ 
nant pd| : 


det{H) < (43) 

Similarly, for the matrix permanent [H 

per(iJ) > per(ij(i’^))per(i?( 2 . 2 ))_ ( 44 ) 

By repeated application of Eqs. (03]) and (El) it is easy 
to demonstrate that Eq. m holds. 


Ill Factorization of J-function and its represen¬ 
tation through a density matrix 


To show that the probability p)^ ') for the com¬ 

pletely indistinguishable bosons (fermions) corresponds 
to the absolute maximum (minimum) over arbitrary in¬ 
put states of particles in a given configuration n, one 
has to know to what class of functions the physical J- 
functions, i.e., describing an input of a linear network, 
belong. Let us show that any normalized by J{I) = 1 
p.s.d. function J(cr) can be represented in the form of 
Eq. (1551) with some state ([551) . Since sgn(CT) J(cr) 
is also a normalized p.s.d function, it is sufficient to con¬ 
sider bosons, e[a) = 1. 

Consider a linear subspace C C defined as the 

linear span of vectors |cr) = Pa\I) for cr G 5jv, where 
|/) = ® ® |(/) 7 v) for some arbitrary orthonormal 

vectors \cj)k) G "H, {(j)k\(t>i) = h,i (note that (ctIt) = 6cr,T 
and Ptt = when restricted to C). Using that 

Pair) = |fTr) and {Tr\Pa\n) = S^j, by starting from a 
trivial identity we get 


M H 


Po\t^) 


= Tr 


(45) 


where the trace is taken in and we have introduced 
a p.s.d. Hermitian operator (density matrix) in the 

Hilbert space 

^ ^ E E (46) 

t^Sn ttGSn 


Obviously, Tr{p(“‘)} = J{I) = 1. Positivity of 
would follow from the explicit form (now for bosons and 
fermions) 


reSw 

|4>r) = ^{crT)P„{\(j)i) . .<^\(j)N)}, (47) 

ctGSn 

where ^{a) = e{a)9*{a) for the factorizing function ff(a) 
(see Eq. (15T1) below). 


Let us sketch the prove of factorization of a p.s.d. Her¬ 
mitian function J{a) (see Ref. [13 )• Consider an opera¬ 
tor J” in £ 

(48) 

<7£Sn 

given by following matrix 

J^,r = {v\J\t) = J{v~^t), (49) 

which, by assumption that J{<j) is a p.s.d. Hermitian 
function, is a p.s.d. Hermitian matrix. Since operators 
having the form given by Eq. (1551) constitute an sub¬ 
algebra of operators in £, the p.s.d. Hermitian operator 
J (l48l) can be factorized by an operator B £ C 

(7 

By the group rule P^Pt = Par Eq. (1501) . in the matrix 
form, is equivalent to the factorization 

'4(cr)= E E = (51) 

t^Sn <t^Sn 

The symmetry (1341) means that the factorizing function 
can be chosen to satisfy 

0{aTr) = 9{a), Vtt G Sn- (52) 

(In contrast, d(Tcr) for all r G is the same factoriza¬ 
tion with a different order of terms). 


IV Derivation of the average probability formu¬ 
lae 

Though the average quantum detection probabilities 
follow from a simple symmetry argument for a Haar- 
random unitary 17, the classical case requires a bit more 
of insight. Here these results are derived by direct eval¬ 
uation demonstrating also the validity of the classical 
formula for M ^ 1 (arbitrary 1C and N), as observed 
in numerical simulations. The following identity will be 
employed 

N 

N 

= E (53) 

a—l 

where yV(M, a) is the Weingarten function of the unitary 
group |14l . [l 6 | which depends only on the cycle structure 
of the relative permutation a = , i.e., the sequence of 

numbers (ci(tT),..., cn{(t)) of cycles of lengths (1,..., N) 
in its cycle decomposition (for more details see Ref. [13 )■ 
By application of Eq. dMl) we have 

N 

H E (54) 

ot—l h'e5n re5m 
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with summation over permutation invariance subgroups 
Sn and iSm of the input and output indices, respec¬ 
tively. Then, from Eq. ra we have for the average over 
Haar-random network (here the output indices vary be¬ 
tween 1 and 1C and, since we sum over the output modes 
li,... ,1 m, one of permutations in J is redundant, giving 
a A^! factor) 


MJ)) = 


N 


M(n)^ M(m) 




(J^Sn 




= E 


' N\ 

^(m) 


E yV{M,aT), 




(55) 


TGSri 


where is over all occupations m of /C output modes 
and we have used the property of J{a) in Eq. (IMl) and 
that the number of G 5n is |5n| = /.t(n). Let us now 
consider separately bosons, fermions and distinguishable 
particles for a general input /i(n) > 1. 

In case of the completely indistinguishable bosons, 
J(ct) = 1, we obtain from Eqs. (l55|) and (l54l) 


m 


^(m) 


E E ^(M,aT) 

(TGSn T^Sm 


{ic + N-iy. 
(/C-l)! 


E 

(t^Sn 


(56) 


where we have used that Sm C Sn and |5m| = (the 
first sum gives the number of Fock states of N bosons in 
/C modes scaled by A^!). Observing that for /C = M the 
probability must be equal to 1, we get the sum of W- 
functions in Eq. (|56ll 


E >V(M,a) 

(tGSn 


(M- 1)! 
(M-h A^- 1)!’ 


(57) 


Combinig Eqs. (IMl) and (E3 we get the final expression 
for 

In case of the completely indistinguishable fermions 
J((t) = sgn((T) and /r(n) = /i(m) = 1 (no two or more 
particles per mode). We have from Eqs. (1551) and (15^ 


(p^^)=Af!E 

(t^Sn 

= TF^W E sgn(a)W(M,a) (58) 

(jGSn 


In case of distinguishable particles from Eqs. (l3^ and 
(1551) we get 

(Pw) = E'--E E 

ll — l /a 7 —1 rG5m 

/c k: N 

= E E E ■ ■ ■ E n 

ly^SnO-^SN h—l /iV —ICK —1 

= E E W{M,v(j)K*'^ (60) 

(T^Sn 

where we have set #cr = Ci (cr) -I- ... -I- cm{(t) (the total 
number of cycles in the cycle decomposition of cr) and 
used that |5„| = p(n) and Oa^i 

Eq. (15(11) must coincide with Eq. (1551) for all particles 
in the same input mode, /i(n) = N\, since the limit of 
distinguishable particles Eq. (l35l) is obtained by making 
identical particles from different modes distinguishable, 
while particles from the same input mode behave as the 
completely indistinguishable bosons. On the other hand, 
for a single-particle input it has a form quite different 
from that of Eq. (1551) for a similar input of the completely 
indistinguishable bosons. Apparently, the general expres¬ 
sion would have quite a cumbersome form. Consider the 
special case of single particles. Using an asymptotic form 
of W for M > I [ill 


yV{M,a) = 


(-1)^ 
(25-2)! 


n(-M9,r.<”)(i + o(T^ 


s!(s-l)!’ 

we obtain from Eq. dini) the leading term in the average 
classical probability as a cycle sum 


<7€Sn s—1 

with (ts = —KMgs) 

1 


N 


s ne 


■Xa) 


(62) 

(63) 


(7^Sn ^ 1 


The cycle sum is evaluated by the generating function 
method (see, for instance. Ref. [13) which satisfies the 
following identity 


T{x) = E Znx^ 

N>1 


= exp 



(64) 


(the first sum is the number of Fock states of N fermions In our case (after evaluation of a table sum [H]) we get 

in K. modes). Setting 1C = M in Eq. (|58)) we obtain lF{x) to be 


E sgn(g)W(M,cr) = ■ (59) 

(J^Sn 

J-B-n ( F'^ 

Eqs. (155)) and ((55|) result in the final expression for (p)^ '). 


.rw = exp|-«/gfc^x 



KM 


(65) 
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For K.M 1, the leading order reads 


^ 1 d^JF(o) ^ (-/CM)^ 

^ “ M dx^ iv! 

due to the fact that the derivative of J-{x), 

dj:{x) -2/CM 

dx VI - 4a;(l + Vl - 4a;) 


( 66 ) 


(67) 


has as a factor a slowly varying function of 
\x\ <C 1 (in comparison with jF(x)). Hence 

= (-/CM)^ [1 - iV(Af - l)/2/CM + ...]. Substi¬ 
tuting Eq. (l66)) into Eq. (l62ll we obtain up to 

a factor {l + 0{^)). 


V Classical computations for the assessment 
protocol are only polynomial in N 


The assessment protocol for the standard version of 
BS is dependent on = per(M), where H is a p.s.d. 
Hermitian matrix built from a submatrix of a network 
matrix U according to Eq. (6) of the main text. Such a 
matrix H can be rewritten as 

M 

i7a./3 = <5a./3 - ^ = (68) 

l=IC+l 

where ki ,are input indices. Below we show that, 
in the dilute limit, NL/M —>■ 0 as iV —>■ oo, with a high 
probability over a choice oi L = M — K. output modes, 
only polynomial in N computations are required to esti¬ 
mate per(i7) of Eq. (l68l) to an arbitrary error e polyno¬ 
mial in N~^ [2^. Since we can select L output modes 
completely arbitrary, one can always find a suitable set 
in a polynomial number of trials. Moreover, an approx¬ 
imation formula to the required polynomial accuracy is 
given (Eq. d75|) below). 

We will use a general formula for the matrix permanent 
of a sum of two matrices M , in our case 

N 

per{H) = l + '^{-iy per($[a,a]), (69) 

r — 1 Cti<...<Ci6r 


where a] is an r-dimensional submatrix of $ built on 
the rows and columns a = (oi,..., a^). Now, since for 
NL/M <C 1 the result on the r.h.s. of Eq. (15^ is close 
to 1 (see Eq. (12) of the main text), we expect that only 
a few terms in this expression are required for estimating 
per(i7) to a polynomial in N~^ error. 

Indeed, let us consider how many terms on the r.h.s. of 
Eq. (I69II are required on average over the Haar-random 
networks. By application of the general result given by 
Eq. (12) in the main text (see also the previous section) 
we have 


(per($[a,a])) 


L(£ + l)-...-(£ + r-l) 

M{M + l)-...-{M + r-l)' ^ ^ 


hence, for r <C VN and r <C '/L we obtain 

{Tr)= Y = (1 + ^)^ (^) > 

(71) 

where e ^ r^[l/7V -|- 1/L] <C 1. 

To estimate per(i£) to an error e = for k > 0, 

we can truncate the sum on the r.h.s. of Eq. HMD at an 
order s satisfying (Tg) <C e. To have a single scale N~^, 
let us concentrate on the case given by M = for 

b = 0(1) and (5 > 0 and L = aN with a = 0(1). In this 
case, averaging over Haar-random U gives (per(i7)) = 
1 — 0{N~^). The truncation order s must satisfy 

i.e., for —>• 00 one must choose s > k/6. Let us set 
s = (k + I)/5 (note that s is an independent of N num¬ 
ber) then for ^ 1 the l.h.s. of Eq. (l72)) is smaller 
by a factor 0{N~^) than the r.h.s.. Now, let us show 
that the number of required computations (flops) Cm in 
the summation of the r.h.s. of Eq. (l6^ to the above 
truncation order s = 0(1) scales only polynomially in N 
as A^ —>■ 00. Since computation of a matrix permanent 
of dimension r requires r2’' operations (flops) by Ryser’s 
algorithm , we can estimate the total number of flops 
as follows. First we note that Eq. (15^ truncated at an 
order s has an equivalent simpler expression (to the same 
accuracy) 

per(i7)Ril-k Y {per($[a,a]) - 1} (73) 

ai<...<as 

(we have taken into account that indices ai < ... < ar 
with r < s are contained in one of ai < ... < as and used 
Eq. (IMl) to sum the terms for each subset oi < ... < a^). 
From Eq. dZSl) we obtain 

where we have used s = (k + l)/(5 for the approximation 
error e = 0{N~'^). 

The estimates in Eqs. dZll) and (1731) (and that of Eq. 
(ITT]) ') have a high probability, since by Chebyshev’s in¬ 
equality (for A^ 1) 

Pr{Ts >e)<^= 0{N-y < 1, (75) 

where we have used Eq. ^ [ii]. Since we select 
L = M — K. output modes for the protocol completely 
arbitrary from the output modes of a network, if it turns 
out that for our choice of L modes for a particular net¬ 
work U we need more than Cn flops, we can select again 
(e.g., by using the Gaussian approximation [2^, the prob¬ 
ability that no suitable choice can be found among 0{N) 
non-intersecting subsets of size L = 0{N) from the to¬ 
tal of M = modes decreases to zero at least as 

0{N-^)). 
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One important example is the test against the distin¬ 
guishable particles, for which one can set k = 2(1 -1- 5). 
Indeed, the difference ~Pn^) between the full quan¬ 
tum and the classical cases (see the main text) reads in 
our case 

{pT - P^N^) = O (1^) = o (lV-1-2^), (76) 

therefore, to distinguish the two cases requires using only 
polynomial computations by the above analysis. 


VI Eflfect of the distinguishability error on prob¬ 
ability pn and its Haar-average 


Consider an input state consisting of single bosons 
from N sources (we will consider uncorrelated sources, 
the case of classically correlated ones is a trivial exten¬ 
sion). The corresponding internal state reads 


=p(i) 0 ...0pW, 


(77) 


where is an internal state of the ath boson. The ideal 
bosonic input corresponds (in the case of only classically 
correlated inputs) to all = \4'){4>\- Hence, let us 
assume that each source produces a boson in a state close 
to a pure state 


= (78) 


with the indistinguishability fidelity defined as = 
{(j)\p^°‘'>\(j)). To derive the scaling law for small variations 
1 < 1 we will assume that sources 

emit bosons with a mean fidelity F. Keeping only the 
first order term in 1 — K we obtain from Eq. GZ]) 

p(“‘) « (|</))(<(<|)®^ (79) 

N 

ct—l 


Taking into account that Tr{(5p^“^} = 0, we obtain from 
the definition (155)) and Eq. (1751) up to the first oder in 
1 - F 


J{a) = Tr{p(“‘)p,,} 1 - (1 - F)[N - ci(a)], (80) 

where ci(tT) is the number of 1-cycles (fixed points) in 
the cycle decomposition of permutation cr (ci(tT) terms 
in the sum in Eq. dZSl) give zero). 

Let us first derive the scaling law for the standard ver¬ 
sion of BS. Substitution of Eq. (1551) into Eq. (1551) gives 
(with = per(7J)) 


Pn{J)-P^^'^-{1-F) 


N 


^P^N^ - E Cl(u) 


7€Sn 


a—1 


( 81 ) 


The second term in the brackets on the r.h.s. of Eq. dUD 
is different from the plain matrix permanent by a fac¬ 
tor at each term, counting how many diagonal elements 
of matrix FI it contains. It can be represented in a com¬ 
putable form by observing that the same counting is done 
by multiplying diagonal elements of by a dummy vari¬ 
able X and application of a derivative w.r.t. a: at cc = 1. 
Thus we have derived 


Pn{J) 


with 


, (-8) 
‘Pn 


-(1-F) 




- ^per{i7(x)}a;=i 


(82) 


H{x) = H + {x- l)diag(iJii,..., Hnn) (83) 

The derivative can be evaluated from the Lagrange repre¬ 
sentation of the polynomial per{iJ(a;)}, thus it requires 
approximating V -|- 1 matrix permanents of the p.s.d. 
Hermitian matrix H(x) at distinct values of x to some 
small error e, which can be chosen inversely polynomial in 
N. In the previous section it is shown that this requires 
only polynomial in N computations. 

Now, let us derive a variant of the scaling law which 
is specifically tailored for the scattershot BS. Assuming 
that M ^ (so that one can use the Gaussian approx¬ 
imation for the elements of the Haar-random unitary ma¬ 
trix U 0) we use the expression (1551) into the average 
probability of detecting all N photons in /C output modes 
(i.e., m = (mi, ..., m/c, 0,..., 0), see also Eqs. (1571) and 
(l36l) ') and recall the formula for the average probability 
(p|^^) (in the dilute limit M ^ N'^): 


N\ 


N 


(wW) = E E 

m ' <7 ^Sn 


m “V J 

(84) 




= (P)l”> - (1 - 


M 


(note that the final result has an appealing form, appar¬ 
ently indicating that it might be valid without assuming 
the dilute limit M ^ V^). Mathematical details in the 
derivation of Eq. (|57)) are as follows. We have used an ap¬ 
proximate independence of 2N matrix elements of U for 
M ^ N'^ and the Gaussian approximation (2^ in com¬ 
puting the average of the product, with {\Uki\^) ~ 1/M 
and for single bosons at input (ka ^ kp for a 13) 


N 








Q!=l 


M^ 




(85) 


TTGSn 


since all permutations of the output indices belonging to 
the output subgroup Sm equally contribute to the result. 
The double sum in the line before the last in Eq. (l84l) 
(over all m and the subgroup 8^) 

« = !«-■=.Ml 

m ' treSm 


( 86 ) 
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is computed as follows. We have a sum over permutations 


1 

/r(m) 






1 

/i(m) 


II H Ilci(az) 

CTjC^^rrif^ 1 


(87) 


i.e., the average number of fixed points (over all permu¬ 
tations) in each permutation group Smi is equal to 1 [la, 
whereas if m; = 0 there is no corresponding contribution. 
The hrst sum in Eq. (I5S1) is then easily computed 


R = 


/ ^ 

hv-/c + ^<5^,,o 


/=1 


= (iV-/C) 
= {N-1) 


{K + N-iy. ^(/C + iV-2)! 
' K.- 


{k: + n-2)i 


iV!(/C-2)! 


(iV- 1)!(/C- 1)!’ 


( 88 ) 


where the summation gives the number of Fock states of 
N bosonic particles distributed over, respectively, K, and 
K. — 1 modes. 


VII Equivalent description of a lossy linear net¬ 
work 

Below we focus on bosonic particles, for fermions 
one should replace the commutators below by the anti¬ 
commutators. A realistic network U is non-unitary due 
to (generally path-dependent) losses of particles, its ac¬ 
tion is described not just by input and output mode op¬ 
erators, oi,..., um and 6i,..., 6 m, but also by some ad¬ 
ditional operators /i,..., /m accounting for losses: 


Eq. (HOI) can be satisfied if we expand the loss operators 
fl,..., over some additional creation operators 

M 

fk — H (91) 

1=1 

where V is a matrix satisfying the following matrix equa¬ 
tion (valid both for bosons and fermions) 

VV'<=I-UU\ (92) 

where / is the unit matrix and denotes the Hermitian 
conjugate to matrix U. Notice that Eq. (ESI) requires 
that the singular values of U be bounded by 1, which 
is the necessary and sufficient condition for an arbitrary 
complex matrix U to describe a passive linear quantum 
network. There is a polar decomposition, U = 
where 77 is a unitary matrix and A = UU^ a Hermitian 
matrix (describing losses in the network) with the eigen¬ 
values bounded by 1. 

The expansion in Eq. dSU means that one can imbed 
an arbitrary (non-unitary) linear M-mode network into 
a 2M-mode unitary one [^. The following embedding 
unitary network seems to be the simplest one 


U = 


U V \ 

-Vtw D J ’ 


(93) 


here the diagonal matrix D = diag(rji,... ,7 ]m), 
0 < r/k < 1, is composed of the square-roots of sin¬ 
gular values of U (eigenvalues of \/A, since A = 
UU^), V = SQ, with the unitary matrix S con¬ 
taining the_eigenvectors_of a/A, VA = SDSy Q = 
diag(-\/l — rjf ,..., and U is from the polar 

decomposition U = -x/AW. Matrix U can be also rewrit¬ 
ten in a product form 


a 


t 

k 


M 




(89) 



D 

-Q D ) 


S'! 0 \ / 77 0 \ 

0 7 y 0 / y ’ 


(94) 


where operators fk and fl commute with the creation 
and annihilation operators corresponding to network 
modes: [fk,ai] = [fk,bi] = [fl,ai] = [fl,bi] = 0 [2l|. 
Using the latter we obtain 

M 

[fkjj]=0, [fkj]] = Sk,j-J2uliU,,i. (90) 

1=1 


from which it is evident that U is unitary. Note that in 
description of a lossy network U, the unitary matrix U 
has vacuum input in the modes {M + 1,..., 2M} and 
output modes {M + 1,..., 2M} are not accessible “loss 
channels”. The embedding matrix of Eq. (1^ reduces 
to a matrix appeared in Ref. in the special case of 
path-independent losses, i.e., a diagonal loss matrix A. 
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